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OO Fluid membranes made out of lipid bilayers are the fundamental separation structure 

Cn in eukaryotic cells. Many physiological processes rely on dramatic shape and topological 

^^ changes (e.g. fusion, fission) of fiuid membrane systems. Fluidity is key to the versatility 

,^ and constant reorganization of lipid bilayers. Here, we study the role of the membrane 

h^- intrinsic viscosity, arising from the friction of the lipid molecules as they rearrange to ac- 

^^ commodate shape changes, in the dynamics of morphological changes of fiuid vesicles. In 

F^ particular, we analyze the competition between the membrane viscosity and the viscos- 

^ ity of the bulk fiuid surrounding the vesicle as the dominant dissipative mechanism. We 

Ph consider the relaxation dynamics of fiuid vesicles put in an out-of-equilibrium state, but 

i_i conclusions can be drawn regarding the kinetics or power consumption in regulated shape 

changes in the cell. On the basis of numerical calculations, we find that the dynamics 

arising from the membrane viscosity are qualitatively different from the dynamics arising 

from the bulk viscosity. When these two dissipation mechanisms are put in competition, 

(T^ we find that for small vesicles the membrane dissipation dominates, with a relaxation 

On time that scales as the size of the vesicle to the power 2. For large vesicles, the bulk 

"^" dissipation dominates, and the exponent in the relaxation time vs. size relation is 3. 
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1. Introduction 

Amphiphilic membranes are self- assembled structures made out of lipids or other am- 
phiphilic molecules such as diblock co-polymers (Discher et al. 1999). Above a transition 
C^ temperature, these membranes are fiuid within the membrane surface, while they retain 

the transversal order, which confers them with bending rigidity. The membrane fiuid- 
ity is essential for many biochemical processes involving membrane proteins (Saffman 
& DelbrucK ly/o). The fiuidity of amphiphilic membranes has been studied experimen- 
tally, with measurements of its two-dimensional viscosity (Danov et al. 2000), as well 
as computationally with molecular dynamics (MD) simulations of viscometric tests (den 
wLier & Shkulipa 2007). Atomistic molecular dynamics (MD) simulations have been very 
useful in this and other respects, but remain limited to small membrane patches due to 
the prohibitive number of atoms involved in closed vesicles and the slow equilibration 
times of membrane systems (v^uuKe ei at. zuuo). Coarse-grained MD allows us to reach 
larger systems, and there has been notable successful studies in recent years involving 
out-of-equilibrium phenomena at the scale of small vesicles, e.g. Reynwar et al. (2007). 
Yet, as in atomistic MD, the computational cost scales as the size of the system to the 
power 6, and therefore there exists a hard upper bound on the sizes of systems that 
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can be reached with current computers. On the other end of the spectrum of model- 
ing approaches for fluid membranes, continuum mechanics has proven very effective in 
describing the mechanics of membrane systems and reproduce experimentally observed 
equilibrium shapes (Lipowsky 1995). There has been much less progress in using con- 
tinuum mechanics models to describe the dynamics of fluid membrane systems, and in 
particular to account for the membrane viscous flow, i.e. the rearrangement of lipids 
on the membrane surface driven by external actions or required to accommodate shape 
changes. This scarcity of results contrasts with the importance of out-of-equilibrium and 
dynamical events both in the cell and in synthetic systems (Baumgart et al. 2003; Bacia 
et at. 'l[)[)b). Here, we address the relaxation dynamics of vesicles (closed membranes 
made out of fluid amphiphilic membranes) with continuum mechanics models and sim- 
ulations, with a particular emphasis on the role of the membrane dissipation due to 
the friction between the amphiphiles as they shear to accommodate shape changes, and 
the competition in setting the relaxation dynamics between this membrane flow and the 
flow induced in the ambient fluid by shape changes. Our simulations can treat without 
difficulties large systems and very slow processes, in contrast with MD methods, at the 
expense of the molecular details, which are crucial for some phenomena. 

The continuum mechanics formulation of the coupling between shape dynamics and 
membrane hydrodynamics was first established by :5criven (lyou) (see also Aris (lyoz)). 
The Navier-Stokes equations on a two-dimensional, time-evolving, curved manifold were 
given in an intrinsic manner in the parameter space of the surface, and required heavy 
differential geometry artillery. We note that the inertial terms in the surface Navier- 
Stokes equations are easily dealt with, and the main complications arise from the viscous 
terms; here we are interested in the low Reynolds number regime, and therefore ignore the 
inertial terms. An alternative cartesian formulation in terms of time-dependent projection 
operators was proposed later (Secomb & Skalak 1982; Barthes-Biesel & Scaler 1985). 
None of these formulations have been exercised beyond surface flows on fixed, simple 
shapes, or beyond infinitesimal shape perturbations around simple shapes such as the 
sphere. For finite shape changes, only extremely restricted families of idealized shapes 
have been considered, e.g. Fischer (1994). Maybe due to the inherent complexity of 
the equations governing the coupled shape dynamics and membrane flow, the effect of 
membrane viscosity has often been neglected in continuum studies, and the idea that 
this source of dissipation can be safely ignored in most situations of interest, compared 
for instance to bulk dissipation, has prevailed. Here, we challenge this viewpoint, which 
of course can be well justified in important situations such as the behavior of vesicles 
in shear flow (P*ozrikidis 2001; Danker & Misbah 2007; Biben et al. 2005). But even in 
this case, recent coarse-grained MD simulations suggest that the membrane viscosity can 
play a crucial role (Noguchi & Gompper 2005a, 6). 

We present, to the best of our knowledge, the first vesicle calculations involving finite 
shape changes and incorporating the effect of the membrane flow (as well as the effect 
of the bulk flow). We perform relaxation dynamics simulations of axisymmetric vesicles 
driven by curvature elasticity, dragged by membrane and bulk viscosity, and subject to 
the usual incompressibility and inextensibility constraints of the bulk and the membrane 
fluids respectively. The new geometric and variational formulation of the Stokes equations 
on a curved, time-evolving surface coupled to the ambient fluid flow and the curvature 
elasticity proposed in /^rruyo (^^ ut^oimune (zuub') is the basis of our numerical strategy. 

The mathematical formulation of the relaxation dynamics of axisymmetric vesicles 
is given in Section 2. The general equations are particularized, and strikingly simple 
expressions are developed for the membrane dissipation functional and the local inexten- 
sibility constraint, which make the theory accessible to analytical calculations and simple 
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numerical implementations. The numerical discretization of the theory, with Galerkin fi- 
nite elements and boundary elements based on B-Splines is described in Section 3. The 
numerical results investigating the qualitative behavior of the dynamics engendered by 
different dissipation mechanisms, and exploring the competition between membrane and 
bulk viscosity are reported in Section 4. The conlcusions we can draw from our analysis 
are collected in Section 5. 



2. Formulation of the relaxation dynamics 

This section describes the continuum model for an inextensible fluid membrane with 
curvature elasticity immersed in a viscous, incompressible Newtonian fluid. We ignore 
inertial forces, and model the amphiphilic membrane as a Newtonian two-dimensional 
fluid, in agreement with coarse-grained molecular dynamics simulations (den Otter & 
Shkulipa 2007) and experimental observations (Dano^/ ^+ ^^ "^^^0; Dimova et ah 2006; 
Cicuta et al. 2007). Besides the effect of the membrane viscosity, which is original to the 
present work, the ingredients in the formulation are classical and well-known, although 
we introduce some particularly simple expressions for axisymmetric vesicles. We firstly 
describe the kinematics, which we choose to describe in a Lagrangian manner both in the 
normal and tangential directions. We then describe the dissipative viscous mechanisms 
in terms of the Rayleigh dissipation potentials, the energetic mechanism operative in this 
model (curvature elasticity), and the constraints. We derive the force balance equations 
from a variational principle, from which the relaxation dynamics follow. 

2.1. Kinematics 

We describe parametrically axisymmetric vesicles in terms of the generating curve, i.e. 
the vesicle surface Vt at a given instant t is given by 

x{u,0;t) = {r {u; t) cos 0,r{u;t) sin 0,z{u;t)}, u G [0,1], G [0,27r], 

where 

c(u] t) = {r{u; t), z{u; t)}, u G [0, 1], 

is the parametric description of the generating curve Ct at the instant t. We consider 
closed surfaces with continuous tangents, hence we require 

r(0) = r(l)=0, z{0)=z{l) = 0, (2.1) 

where (•)' denotes partial differentiation with respect to u. 

For simplicity, from this point on we omit the dependence of all quantities on time. We 
shall formulate the mechanics of the fluid membrane in terms of the generating curve. 
Its speed is given by a{u) = y^^PT^^)P^^T^^T^^OF- Integrals on the surface can be brought 
to the interval [0, 1] with the relation d^ = {27Tar)du. The tangent unit vector to the 
generating curve pointing in the u direction and a unit normal are given by 

t=^{r',z'}, n=^{-z',r'}. 

At this point, there are several possibilities in describing the kinematics of the mem- 
brane. Since amphiphilic membranes are two-dimensional fluid bodies moving in a higher- 
dimensional space, the description of the motion is necessarily Lagrangian, at least with 
regards to the shape of the membrane, i.e. c-n = Vn onT, where the dot denotes partial 
differentiation with respect to time and Vn is the normal velocity of the membrane. In 
the tangential direction, however, there is gauge freedom due to the particle relabeling 
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symmetry of the fluid flow equations, here the re-parameterization invariance with re- 
spect to u. One special choice is the Eulerian gauge, for which c • t = and Vt, the 
tangential velocity of the material particles on the membrane, enters the formulation as 
an independent variable. Here we adopt another natural choice, the Lagrangian gauge, 
in which the parameter u is viewed as a label for material particles, hence c - 1 = Vf. In 
this case, we have: 

C = {r^ z} = V = V -\- VnTl = Vtt + VnTl^ 

and 

vt = -{r'r + z'z), Vn = -{-z'r + r'z). (2.2) 

a a 

This decomposition is important since only the normal velocities change the shape of the 

membrane, while the tangential velocities represent the flow of the amphiphiles on the 

membrane surface. The above relations can be inverted as 

r = -{r'vt - z'vn), z = -{r'vn + z'vt). 

A key object for the kinematics and constitutive relations of the two-dimensional fluid 
on a curved, time-evolving surface is the rate-of-deformation tensor d. This tensor can 
be geometrically thought of as the tangent projection of the rate of change of the metric 
tensor of the surface as it is advected by the membrane velocity (Scriven 1960; Marsden 
& Hughes 1983; Arroyo & DeSimone 2009), which leads to the expression 

d = - {VsV + VsV^) - Vnk, 

where Vg denotes the surface covariant derivative and k = — Vgn the second fundamen- 
tal form. The formulas for axisymmetric surfaces can be found in Arroyo & DeSimone 
(2009). 

Introducing b{u) = —r"{u)z'{u) + r'{u)z"{u)^ we can write the mean and gaussian 
curvatures as 

H=lfl + ^), A- = ^^ (2^3) 

a ya^^ r J a^r 

Note that here H is defined as the trace of the second fundamental form. 

2.2. Governing equations 

We present first the generic ingredients entering the governing equations for the shape 
evolution, under the assumption of low Reynolds number, hence neglecting inertia. The 
Rayleigh dissipation potentials can be expressed by a functional W[r^ i], which, although 
not explicitly written, depends nonlinearly on the shape of the surface. Under the as- 
sumption that both the ambient bulk fluid and the membrane two-dimensional fluid are 
Newtonian, this functional is quadratic in its arguments. 

The curvature energy and in general other energetic mechanisms such as line tension, 
depend exclusively on the shape of the membrane, hence can be expressed in terms of 
a nonlinear functional H[r, z]. The energy release rate functional is the negative of the 
variation of the energy functional in the direction of {r^z}^ i.e. G[f^z\ = — (5H[r, 2:; r, i], 
where again its explicit (nonlinear) dependence on the configuration of the surface has 
been omitted. By its definition, it is clear that the energy release rate functional is linear 
in its arguments. 

The membrane dynamics are often constrained by global nonlinear equalities such as 
enclosed volume or surface area constraints. We generically express these constraints with 
a vector-valued functional as B[r^z] = 0, or in rate form as C[r^z\ = 5B[r^z]f^z\ = 0. 



The role of membrane viscosity in the dynamics of fluid membranes 5 

Again, this functional is linear in its arguments and nonlinear in the configuration of the 
membrane. We shall also consider local constraints such as the local inextensibility of 
the membrane, expressed at each point of the membrane as c{r^ z) = 0. 

At each instant, the dynamics of the membrane equilibrate the dissipative and the 
energetic forces subject to the constraints. In other words, the membrane evolves in 
such a way that the rate of viscous dissipation resisting the fiow is exactly equal to the 
energy release rate during the fiow. Mathematically, the evolution equations follow from 
minimizing iy[r,i] — G[r^z] subject to the constraints. Forming the Lagrangian 

£[r, i. A, A] = W[r, z] - G[r, z] - f Xc{r, z)dS - A • C[r, i], 

the velocities and the Lagrange multipliers at each configuration {r, z} are obtained as 
the stationary points with respect to all admissible variations 

5rC = 5iC = 5xC = 5aC = 0, (2.4) 

which is a form of the principle of virtual power. Note that, unlike the other arguments 
of the Lagrangian functional, the Lagrange multipliers for the global constraints A are 
not functions of u. If c{r^ i) = expresses the local inextensibility of the membrane, then 
the Lagrange multiplier A is the surface tension. We provide below the specific form of 
these functionals. 

2.3. Constraints 

Fluid membranes are semipermeable, and assuming osmotic equilibrium between the 
enclosed and the outer media, it is often reasonable to assume that the enclosed volume 
is constant (Seifert & Lipowsky 1995). This constraint is expressed as 

= C^°n^. z]=V = - [ VndS = - [ -{-z'r + r' z){2iTar)du. (2.5) 

Jv Jo ^ 

Generally, the in-plane stresses on the membrane are low compared to the elastic moduli, 
and the membrane can be assumed to be locally inextensible (Seifert & Lipowsky 1995). 
This condition on the surface can be expressed as (Arroyo & DeSimone 2009) 

= trace d = Vg • v — Hvn, 

which for axisymmetric surfaces reduces to 

= c'^'^^'ir.z) = ^{rvtY - Hv^. (2.6) 

ar 

With the Lagrangian gauge mentioned earlier, this constraint takes the particularly sim- 
ple expression (see appendix B for a derivation) 



9 , a r 1 



r 



^(r, i) = ^ In ar = - + - = ^{r'r' + z'z') + -. (2.7) 

at a r a^ r 

If the details of the fiuid fiow on the membrane are disregarded, this local condition is 
often replaced by a total surface area constraint (uu ei at. zuuy), which for a closed 
surface takes the form 



= C^'^^[r, z]=S = [ trace d dS = - [ Hv^dS = - [ —{-z'r + r'z){27Tar)du, 
Jv Jv Jo ^ 

(2.8) 

where the divergence theorem has been used. 



6 M. Arroyo, A. DeSimone and L. Heltai 

2.4. Membrane dissipation 

As derived in Arroyo & DeSimone (2009), the Rayleigh dissipation potential for a New- 
tonian closed fluid membrane can be written as 



[V,Vn\ 



ji d : d dS 



L 



/i 



dv 



b|2 



(V, . vf -K\v\^^ {H^ - 2K)vl - 2{VsV : k)v^ 



dS. 



where d denotes the exterior derivative, dv^ is the generalization of the curl of the tangent 
velocity fleld on the surface, and /i denotes the membrane viscosity, with units of force x 
time X length"-^. Molecular dynamics simulations (aen utter iSL :::)nKunpa zuu/) as well as 
experiments (Danov et al. 2000; Dimova et al. 2006) support modeling fluid membranes 
as Newtonian two-dimensional fluids. For axisymmetric surfaces, this expression reduces 
to 
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dS. 



The matrix in the above expression can be checked to be positive semideflnite with rank 
2, which reflects the zero-dissipation mode resulting from translating rigidly the vesicle 
along the symmetry axis in the z direction, as a consequence of the internal nature of 
the membrane dissipation. 

We provide here a remarkably simple alternative expression of this dissipation potential 
under the hypothesis of axisymmetry with the Lagrangian gauge. We flrst use the the local 
inextensibility constraint to express v'^ in terms of Vt and Vn- The dissipation potential 
then takes the form 



W 



I 



1^-, 



[-* 



" {ar)'^ 
which, using equation (2.2), reduces to 



{r') 



i\2 



[z'Y 



Vt 
Vn 



dS. 



2/i 



dS 



2/i 



(27rar)dii. 



(2.9) 



2.5. Bulk dissipation 

We assume that the velocities of the membrane and of the surrounding fluid coincide, 
i.e. no slip as suggested by coarse-grained MD simulation (den Otter & Shkulipa 2007) 
and conventionally assumed for lipid- water interactions (Stone & Ajdari 1998). This 
hypothesis may not be adequate in extreme situations, for instance if there is signiflcant 
flow across the membrane. We consider an inflnite incompressible fluid at rest at inflnity 
and axisymmetry. 

In the absence of body forces, the Rayleigh potential accounting for the dissipation 
in the surrounding fluid can be transformed into a surface integral by virtue of the 
divergence theorem as 



W 



bulk 






DdV 



Ag • VAS, 
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where a is the stress tensor and D the rate-of-deformation tensor in the bulk fluid, V de- 
notes the velocity field on the membrane, and Ag denotes the jump of tractions across the 
membrane due to the elastic forces, constraint reactions, and membrane viscous forces. 
By inverting the boundary integral representation of the velocities on the membrane in 
terms of a single layer potential, which assumes that the fiuid viscosity is the same inside 
and outside of the vesicle, we obtain 

^N = -TtL / G(cW, c{v)) . Ag{v) dS{v) := ^(AAg){u), (2.10) 

where G is the axisymmetric Green's function and ji^^^^ is the bulk fiuid viscosity. This 
representation allows us to implicitly express ly^^^^ in terms of r and i. 

The axisymmetric kernel G represents the single layer integral of a ring of singularities 
passing through the point c{u): 

^27r r 

Q{c{u)^c{v)) := r{v) 



/o 
where S is the three-dimensional Stokeslet 

1 / I 



^zz (S^^ COS 6> + Sa,^ sin 6>) 

^zy {^yy COS + S^,^ siu 0) 



di9, (2.11) 



^^'^" 8.VM 

evaluated at r = R{0)c{v) — c{u)^ with R{0) a rotation of angle around the axis of 
symmetry. 

A detailed derivation of equation (2.11) can be found in Pozrikidis (1992), together 
with an explicit expression of G in terms of complete elliptic integrals of the first and 
second kind. Using the implicit definition of the Dirichlet to Neumann map A~^ given 
in equation (2.10), we can then express the bulk dissipation as 

W^^^^[r, z]= [ -/i^^^^ V • A'WdS = [ -/i^^^^jr, i} • A'^ir, z}{27rar)du. (2.12) 
Jr 2 Jo 2 

This expression is valid under the compatibility condition J^V ' n dS = 0, a result of 
the fiuid incompressibility, enforced by Eq. (2.5) 

2.6. Curvature elasticity 

The curvature elasticity of the fiuid membrane is modeled with the Helfrich-Canham 
functional (see Seifert & Lipowsky 1995, for a discussion on curvature elasticity models), 

U= f '^{H - Cof dS^ f ncK dS, (2.13) 

where Co denotes the spontaneous curvature and n and kq are elastic parameters. For 
closed homogeneous surfaces, the second term is a topological invariant, which we will 
disregard here. For a more general treatment including membranes with boundary, see 
Arroyo & DeSimone (2009). The energy release rate takes the form 



G[Vn] = - J k{H - Co) 



A,v„ + ^{H^-4K + HC'o) v„ 



dS. 



For the numerical implementation, to avoid third order derivatives of r and z, it proves 
more convenient to take variations directly from the following form of the elastic energy 

/•I ^ 
n[r,z]= -{H-Cofi2TTar)du, (2.14) 

Jo ^ 
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together with equation (2.3), which yields an expression for G[r^z] given in Appendix 
A. Both approaches can be shown to be equivalent by integration by parts with the 
conditions in equation (2.1) to annihilate the boundary terms. 



2.7. Different models 

In the paper, we consider several models to assess the specific effects of the membrane 
viscous flow and its relevance. Initially, for the sake of qualitative comparison, we consider 
a model with only membrane dissipation, and local membrane inextensibility 

£^[r, i. A, A^°^] = W"^^"^[r, z] - G[r, z] - f Ac^"^^*(r, z)dS - A^°^C^°^[r, i]. 

We consider also a model with only bulk dissipation, together with a global area con- 
straint as in uu ei at. (zuuy) 

C^[r, i, A^^^% A^°^] = W^^^^[r, z] - G[r, z] - A^^ea^areaj^^ ^] _ Avol^volj^^ -^^ 

A similar model with local inextensibility can be considered, as in Biben et al. (2005). 
We have checked that the results for both models are very close, hence consider only C^ 
for definiteness. 

We also consider a model whose dissipative mechanism is of purely mathematical intent 

C^[r, i, A^^^% A^°^] = W^' [r, z] - G[r, z] - A^rea^areaj^^ -^ _ A-oiC"°i[r, i], 

where 

W'^^[r,z]=l^^vldS, 

and /i is a mathematical viscosity coefficient. The resulting dynamics, for the closed 
surfaces considered here, are a constrained version of the L2 gradient flow of the Willmore 
energy. Such a gradient flow finds applications in geometrical analysis (Simonett 2001), 
and is also considered by some authors as a simple model for the dynamics of fluid 
membranes {uu a at. zuuo). In all the models presented this far, vesicles of different 
sizes evolve in the exact same way, upon re-scaling of the time variable. Finally, the 
following model accounts for both the membrane and the bulk dissipation 

C^^'^^r, i. A, A"°^] = V^"^^"^[r, z] + W^^^^[r, z] -G[r,z]- f Ad^"^"*(r, z)dS - A"°^L)^°H^, ^]- 

In this model, the two dissipative mechanisms compete, and vesicles of different sizes 
exhibit a different relaxation behavior, as reported in ^nu^u 06 DeSimone (2009) and 
illustrated later in the present paper. 

These models could be supplemented by external actions, although here we only con- 
sider the relaxation dynamics from an initial out-of-equilibrium condition. The out-of- 
equilibrium configurations are obtained as equilibria for a given value of the spontaneous 
curvature Co, which in the dynamics simulation is set to zero, or by applying an exter- 
nal force on the vesicle, which is suddenly released so that the systems returns to an 
equilibrium state. It should be mentioned that, while for models £^, £^, and C^^^^ the 
boundary conditions in equation (2.1) are sufficient, the model C^ requires an additional 
condition to fix the motion along z, which otherwise remains undetermined due to the 
invariance of 1^"^®"^ apparent from equation (2.9). 
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3. Numerical approximation 

3.1. Spatial semi- discretization 

The spatial discretization follows from a standard Galerkin approach. The generating 
curve of the axisymmetric surface is represented numerically as a B-Spline curve 

N 

c{u;t) = {r{u;t),z{u;t)} ^ V5/(ix) {r/(t), z/(t)}, 

— ^ .. ^ 

where Bi{u) are the B-Spline basis functions (Piegl & Tiller 1997) defined on the interval 
[0, 1], and {r/(t), z/(t)} is the position of the /— th control point of the B-Spline curve at 
instant t. Again, we drop the dependence on time. The velocity of the membrane can be 
computed as 

N 

V(«)«^i?,HP,. (3.1) 

1=1 

Since the formulation involves up to second derivatives of the generating curve, it is 
convenient that the numerical representation is sufficiently smooth, hence avoiding cum- 
bersome mixed approaches. Here, we have considered cubic B- Splines, which have up to 
second-order continuous derivatives. With the natural numbering of the basis functions, 
the symmetry conditions in equation (2.1) can be expressed in this numerical representa- 
tion as ri = rjy = 0, zi = Z2 and zn-i = zjy. We collect all the nodal values in the vector 
P = {ri, ^1, r2, 2:2, . . . , Tat, ^at}. If local inextensibility is required, we need to discretize 
the field of Lagrange multipliers 

M 



X{u)^Y.BjXj, 



j=i 

where Bj are taken here to be the quadratic B- Splines obtained from the same knot 
vector used to generate the functions Bj. We have checked numerically the stability and 
convergence of this velocity /pressure pair in the present context. 

Plugging these representations into the different models, generically written as £, 
and equating to zero the derivatives of jC with respect to velocity of the control points 
{^/(t), i/(t)}, possibly Xj and the global multipliers A (or equivalently evaluating the 
variations in equation (2.4) at the B-Spline basis functions), the Galerkin semi-discrete 
equations follow as 

D(P)P + L(P)A = /(P) (3.2) 

L^(P)P = 

where A collects all the Lagrange multipliers for a given model. 

Specifically, recalling equation (2.5), the entries for the column of the constraint matrix 
L corresponding to the global volume constraint are 

Similarly, recalling equation (2.8), the column of L corresponding to the global area 
constraint is 
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The matrix entries of the local inextensibility constraints follow from Eq. (2.7): 

/-w = - / (^B'l + -Bi) Bj{2Trar)du, l}jTf = - [ ^B'iBj{2Trar)du. 
Jo \^ ^ J Jo ^ 

Recalling equation (A 1), the vector of nodal elastic forces is computed as 

fir = G[Bj,0], fi, = G[0,Bj], 
From equation (2.9), it follows that the only nonzero entries of D"^®"^ are 



DTrTr= [ ^BjBj—du. 

Jo ^ 



3.2. The hulk dissipation matrix D^^^^ 

To discretize the bulk dissipation matrix, we follow a Galerkin Boundary Element method 
based on B-Splines. The traction jumps are discretized using the same B-Spline basis 
functions 

N 

A^(^)^^5,(ii)A^,. (3.3) 

1=1 
Recalling equations (3.1) and (3.3), multiplying equation (2.10) by a test function Bi 
and integrating over the surface, we have 



5ij 



j BiBjdS Pjj = -^ (J Bi{u) J G,k{c{u),c{v))BK{v)dS{v)dS{u)\ AgKk. 



where summation over repeated indices is implied, and the lower-case indices i, j and k 
run over r and z. Let us write the discretized dissipation potential for the bulk fluid. We 
introduce the vector 

G = {Agir, A^i^, A^2r, A^2^, . . . , A^TVr, A^at^} 

collecting the nodal traction jumps. It must be carefully noted that, unlike the vector of 
nodal forces / and the other terms in the first Equation in (3.3), G is not power- conjugate 
to P. We have 



lybulk 

2 



^J^V'AgdS= ^ EE (/^^^^ ^^) ^9i • Pj = Ig^^P. 



which, recalling that G = A ^MP , becomes 

2 
and we finally identify D^^^^ = MA"^M. 

3.3. Implementation details 

The computational strategy described above has been implemented within Matlab, ex- 
cept for the boundary integral method, which has been coded in C++ and is built on top 
of the deal. II library (oan^t^nn ti at. zuu/) with B-Spline support based on the GSL 
library (Gough 2009). A Matlab interface of this code has been set up. This system of 
differential- algebraic equations in (3.3) is integrated in time with the specialized solvers 
odel5s and ode23t. 

Special care must be taken in performing the integration of the diagonal terms of the A 
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Figure 1. Non-uniform knot span adapted to the curvature of the generating curve (left), and 
control polygon for a given vesicle shape (right). The the parameterization of c(u) is nearly 
arc-length. The clustering in the knots of the B-Spline results in a clustering of the control 
points in highly curved regions. 

matrix, since the axisymmetric Stokeslet (2.11) presents a logarithmic singularity when 
its arguments are close to each other. We use integration formulas based on weighted 
Gauss quadrature rules. 

3.4. Reparameterization of the B-Spline curve 

As the system in equation (3.3) is integrated in time, the control points of the generating 
curve collected in P may cluster in some areas, over-resolving some parts of the domain 
and leaving others with poor resolution, irrespective of the geometric features of the 
generating curve. Furthermore, as illustrated in the examples later, the highly curved 
parts of the generating curve may move significantly in the parametric domain [0, 1]. For 
these reasons, it is important to re-parameterize the generating curve periodically in the 
simulations. 

Here, we have devised a heuristic method that builds a non-uniform knot-span for the 
B-Spline functions, which clusters knots, hence basis functions, where the curvature of 
the generating curve is high. At the re-parameterization instants, the new control points 
relative to the new set of basis functions are obtained by a least square fit to the previous 
description of the generating curve that results in a nearly arc-length parameterization. 
See Figure 1 for an illustration. This method allows us to reduce substantially the number 
of control points as compared to resolving the fine geometrical features with uniform 
refinement, and still obtain very accurate solutions. 



4. Results 

Snapshots from a typical simulation are shown in figure 2. The tangential velocity on 
the membrane is represented to visually illustrate the inextensible flow of amphiphiles 
on the membrane, whose tangential shear causes the membrane dissipation. The ambient 
fluid velocity field has also been represented, as post-processed from the boundary integral 
numerical method. This bulk velocity field is incompressible, its tangential projection on 
the membrane is inextensible, and its normal component on the membrane illustrates 
the rate of change of the shape. 

We next report two different sets of relaxation dynamics simulations. On the one hand, 
the qualitative differences of the flows engendered by the different dissipative mechanisms 
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Figure 2. Snapshots of a typical simulation. Sector of the axisymmetric shape (green), tangen- 
tial velocity field on the membrane representing the flow of amphiphiles as they rearrange to 
accommodate the shape evolution (blue), and velocity field of the ambient fluid post-processed 
from the boundary integral solution (magenta). 
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Figure 3. Discocyte vesicle evolving towards the sphere after the volume constraint is released. 
Snapshots of the evolution (top) , and elastic energy of the membrane as a function of time in 
linear and semi- logarithmic scales, for models C"^ and C^ (bottom). 



are analyzed. On the other hand, we consider the competition between the bulk and the 
membrane dissipative mechanisms in C^^^^ as a function of the system size, and obtain a 
composite power law for the relaxation time as a function of size. 

4.1. Comparison of the dynamics for different dissipation mechanisms 

We consider in this section a gallery of relaxation simulations, and compare qualitatively 
the dynamics obtained with each individual dissipation mechanism, i.e. the membrane 
dissipation in model >C^, the bulk dissipation in model £^, and the L2 dissipation in 
model C^ . In the examples in this section, the time scale is arbitrary and has been re- 
scaled so that, for a given test case, all the models reach a given energy threshold at the 
same instant. 

4.1.1. Relaxation of a discocyte vesicle 

We consider first the relaxation of a discocyte vesicle, in equilibrium for zero sponta- 
neous curvature and a reduced volume of 

V = ^"^^^ = 0.65, 

where V is the volume of the vesicle and S its surface area (Seifert & Lipowsky 1995). 
The reduced volume is 1 for a sphere, and smaller for any other shape. In this somewhat 
artificial example, the volume constraint is suddenly released, and the system evolves from 
the discocyte configuration towards the sphere adopting oblate configurations along the 
way. For this example it does not make sense to consider the bulk fluid viscosity, i.e. 
functional £^, since the enclosed volume cannot be incompressible. Figure 3 illustrates 
the relaxation dynamics, as well as the time evolution of the elastic energy. This evolution 
is presented also in semi-logarithmic scale, which better shows the qualitative behavior 
in that it is easy to visually detect deviations from the exponential relaxation of a simple 
linear system characterized by a single time-scale. In this example, model C^ without 
volume constraint behaves nearly exponentially, as shown by the linear response of this 
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Figure 4. Time evolution of the first three eigenvalues of the matrix W for the relaxing 

discocyte vesicle. 

model in the semi- logarithmic plot. The sharp deviation from the exponential behavior 
of model £^ without volume constraint is evident from the plot, particularly at early 
stages. As a matter of fact, the energy evolution for the membrane dissipation model 
seems to have a vertical asymptote at the origin. 

To investigate the origin of this behavior, we eliminate the Lagrange multipliers from 
equation (3.3) by isolating P in the first equation, substituting in the second to isolate 
A, and substituting back in the first equation, which results in 

P = (D-i - D-iL(L^D-iL)-iL^D-i) /(P). 

^" V ' 

W(P) 

The largest eigenvalues of the matrix W {P{t)) characterize the effective inverse viscosity 
of the dominant modes of the system, once the constraints have been eliminated. While for 
the model based on the L2 dissipation, >C^, the eigenvalues remain nearly constant during 
the time evolution, figure 4 shows that the maximum eigenvalue for the model based on 
the membrane dissipation, i2^, changes by three orders of magnitude. This shows the 
strong geometry dependence of the membrane dissipation, which seems to nearly vanish 
for discocytes, and then rapidly increases in the transition to oblate morphologies. We 
revisit later this behavior in the example of a relaxing stomatocyte. 



4.1.2. Relaxation of a pearled vesicle 

We now consider the relaxation of a pearled vesicle with the enclosed volume constraint 
in place. The initial configuration is in equilibrium for a non-zero value of the spontaneous 
curvature. In the dynamical simulation, Co is set to zero, and then the system evolves 
towards a prolate configuration, see Figure 5 (left). The experiment can be viewed as the 
large scale shape evolution after two dissimilar vesicles have fused and the fusion pore has 
expanded to some extent. The energy relaxation with the three models, >C^, i2^, and £^, 
is shown in Figure 5 (right). For this example, in contrast with the example above, the L2 
evolution is the one which deviates most from the exponential relaxation, with an initial 
very fast regime, for which the system can significantly lower its energy without much 
dissipation, then a slowdown and finally an exponential tail. The membrane and the bulk 
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Figure 5. Pearled vesicle evolving towards a prolate configuration, modeling the large scale 
relaxation after a fusion event between two dissimilar vesicles. Snapshots of the evolution (left), 
and elastic energy of the membrane as a function of time in semi-logarithmic scale, for models 
£^, £^, and £^ (right). 



dissipation evolutions deviate slightly from the exponential relaxation, the former with 
a slower evolution and the latter with a faster evolution at early times. 

4.1.3. Relaxation of a tethered vesicle 

We now study the relaxation dynamics of a tethered vesicle, prepared computationally 
following the experimental setup in j_jt-u a tu. (^uu ) in which a vesicle is pulled by 
two beads controlled by optical tweezers. A breaking of the symmetry was observed 
experimentally, and also in our simulations, as in one of the attachments the vesicle 
produces a tether to relieve the elastic energy while maintaining the volume and area 
constraints. Here, we report the relaxation as the forces stretching the vesicle are released, 
see Figure 6. Not suprisingly, the features of the relaxation dynamics are similar to those 
in the previous example. However, in this case, two exponential regimes (linear in the 
semi-logarithmic scale) can be observed, corresponding to the relatively fast event of 
the tether retraction, and to the slower event of the global convergence to the prolate 
configuration. 

4.1.4. Relaxation of a stomatocyte vesicle 

We consider here the relaxation of a stomatocyte vesicle, prepared by minimizing the 
elastic energy with a non-zero spontaneous curvature. The experiment can be viewed 
as the large scale shape evolution after two dissimilar vesicles, one being inside of the 
other, have fused and the fusion pore has expanded to some extent (see Figure 7). This 
last example exhibits a dramatic dependence of the relaxation dynamics on the dissipa- 
tive mechanism. Generically, we observe a very slow dynamics at early stages, where the 
elastic forces are nearly perpendicular to the constraint manifold. An exponential relax- 
ation is observed for the three models. Once the neck of the stomatocyte is erased (third 
snapshot in Figure 7, much faster dynamics are observed, with well-defined exponential 
tails characterized by much larger time constants in both the bulk dissipation and the L2 
dissipation models. In sharp contrast, the membrane dissipation model exhibits an even 
faster, non exponential relaxation at later stages, and abruptly relaxes to the equilibrium 
discocyte morphology. As a matter of fact, there seems to be a vertical asymptote in the 
energy vs time relation, a similar behavior to that observed in the first example, also near 
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Figure 6. Tethered vesicle evolving towards a prolate configuration, after the force keeping the 
vesicle stretched is released. Snapshots of the evolution (top), and elastic energy of the membrane 
as a function of time in semi-logarithmic scale, for models C , £^, and C (bottom). 




?3^ 

?10° 



o^ 10 



-bulk diss 
-L2 diss 
-mem diss 




100 200 300 

time 



400 



Figure 7. Stomatocyte vesicle evolving towards a discocyte configuration, modeling the relax- 
ation after the fusion of two dissimilar vesicles. Snapshots of the evolution (top), and elastic 
energy of the membrane as a function of time in semi- logarithmic scale, for models C , £^, and 
C^ (bottom). 
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a discocyte configuration. Ttiis vanistiing effective membrane viscosity for such configu- 
rations is worth analytical investigation, and highlights the strong geometry dependence 
of the membrane dissipation. 

4.2. Competition between membrane and bulk viscosity 

As a final more realistic example, we consider the relaxation dynamics when both the 
membrane and the bulk dissipative mechanisms are in place. In this case, a natural 
length-scale arises as the ratio between the membrane two-dimensional viscosity and the 
bulk viscosity (Arroyo & DeSimone 2009) 

This length-scale arises also in the Saffman-Delbruck theory for the diffusion of membrane 
inclusions (Saffman & Delbruck 197-^), see also (Stone & Ajdari 199( ), and separates 
the behavior of large vesicles, for which the bulk viscosity dominates over the membrane 
viscosity, and that of small vesicles, whose dynamics are mostly dictated by the membrane 
dissipation. On the basis of a simple analytical model for budding, it has been estimated 
(Arroyo & DeSimone 2009) that for large vesicles, the relaxation time scales as 

..bulk r)3 
tr - ^^^, (4.2) 

while for small vesicles, we have 

*-^- *") 

The cross-over between the two power-laws are for vesicles sizes of Rq ~ ^/2, where Rq 
is the characteristic size of the vesicle, here the radius of the half- sphere whose surface 
area is that of the vesicle under consideration. 

We test these analytical estimates with a simulation resolving the actual shape relax- 
ation coupled with the membrane and the bulk dissipations, as described by the functional 
£fuii ^g perform relaxation experiments of pearled vesicles of different sizes converging 
to an oblate configuration, and for each simulation, we record the relaxation time, de- 
fined as the time it takes for the system to relax up to 99.25% of its excess elastic energy. 
The results are shown in a logarithmic scale in Figure 8, together with the predictions of 
Eqs. (4.2-4.3). The data of the simulations show a remarkable agreement with the com- 
posite power-law predicted by the analytical estimates, not only qualitatively, but also 
quantitatively in this example, even though the budding model in Arroyo & DeSimone 
(2009) differs significantly from the relaxation simulations considered here. 



5. Conclusions 

We have presented, to the best of our knowledge, the first calculations for the dynam- 
ics of fluid membranes accounting for the coupling between the morphological dynamics 
for finite shape changes and the membrane two-dimensional curved flow. To this end, 
we have particularized the theory in Arroyo & DeSimone (2009) to the axisymmetric 
case, yielding a simple, workable model, which is the basis to the numerical simulations. 
The reported simulations show that the relaxation dynamics driven by curvature elas- 
ticity with the usual constraints are qualitatively different for the different dissipation 
mechanisms (membrane viscosity, bulk viscosity and a mathematical L2 dissipation) con- 
sidered separately. Furthermore, when the membrane and the bulk viscosity are operative 
simultaneously, the size of the system dictates the dominant dissipative mechanism. For 
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Figure 8. Composite power law for the relaxation time as a function of vesicle size for a fusion 
event, with a regime dominated by bulk viscosity (large vesicles) and another dominated by 
membrane viscosity. 



vesicles smaller than a given length-scale, membrane viscosity dominates the dynami- 
cal behavior. This length-scale depends on the membrane viscosity, which can exhibit 
variations of several orders or magnitude depending on the amphiphiles chemistry and 
composition, as well as on temperature (Danov et at. 2000). In a significant set of situ- 
ations, the regime dominated by membrane viscosity turns out to be the relevant one, 
as the separation length-scale is on the order of a few to tens of microns for usual lipid 
membranes at physiological temperatures (Dimova et al. 2006), reaching hundreds of 
microns for liquid ordered phases (Bacia et al. 2005) and several millimeters for polymer 
vesicles (Discher et al 1999; Dimova et al 2002; Zhou & Yan 2005). 
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Appendix A. Partial derivatives of the mean curvature 

Taking variations of equation (2.14), it follows that 






r ., z' 



n{H - Co) — r H z + ar irdu, 



(Al) 



where 

Si:^ _ z' dH _ z" 3br' dH _ z' dH _ r" 3bz' 1 dH _ r' 

dr ar'^ ' dr^ a^ a^ ' dr^^ a^ ' dz^ a^ a^ ar ' dz^^ a^ 

Appendix B. Local inextensibility constraint in the Lagrangian gauge 

We rewrite the axisymmetric form of the inextensibility constraint in Eq. (2.6) under 
the Lagrangian gauge given by Eq. (2.2): 
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